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We propose a many-particle- inspired theory for granular outflows from a hopper and for the escape 
dynamics through a bottleneck based on a continuity equation in polar coordinates. If the inflow 
is below the maximum outflow, we find an asymptotic stationary solution. If the inflow is above 
this value, we observe queue formation, which can be described by a shock wave equation. We also 
address the experimental observation of intermittent outflows, taking into account the lack of space 
in the merging zone by a minimum function and coordination problems by a stochastic variable. 
This results in avalanches of different sizes even if friction, force networks, inelastic collapse, or 
delay-induced stop-and-go waves are not assumed. Our intermittent flows result from a random 
alternation between particle propagation and gap propagation. Erratic flows in congested merging 
zones of vehicle traffic may be explained in a similar way. 



Driven granular media display a rich spectrum of pat- 
tern formation phenomena. This includes collective oscil- 
lating states, convection patterns, the spontaneous segre- 
gation of different granular materials, and the formation 
of avalanches due to self-organized criticality Here, 
we will focus on jamming and clogging phenomena Q 
related to arching and intermittent outflows through 
hoppers 0, Q . Similar phenomena are known from dense 
pedestrian crowds |(J. The escape dynamics of individ- 
uals from a room has been intensively studied, showing 
that in crowd stampedes, rooms are emptied in an irreg- 
ular, strongly intermittent fashion ^lj. This effect has 
been discovered in simulations performed with the social 
and the centrifugal force model I 1 11 1 121 , with cellular au- 
tomata and lattice gas automata [13( , and in a meanfield 
model |l4j . It has also been experimentally confirmed 
0, Q- However, analytical models of escape dynamics 
and granular bottleneck flows are lacking. 

In this Letter we will formulate such a model. Our goal 
is to gain a better understanding of (i) the resulting den- 
sity profiles and (ii) the irregular outflows at bottlenecks. 
This includes not only the distribution of the avalanche 
sizes in the outflow from a bottleneck. We will also offer 
a possible explanation of the long-standing problem of 
perturbations forming in merging zones of freeway traf- 
fic flows 0,0], which are characterized by erratic, for- 
ward or backward moving shock waves |l5| . It is believed 
that these can trigger stop-and-go waves in traffic flows 
|l5lll7j . Similar findings have been made in over-crowded 
pedestrian flows J\ and expected for merging flows in ur- 
ban traffic and production networks. 

In all these cases, the competition of too many entities 
for little space leads to coordination problems. We are 
therefore looking for a minimal, common model captur- 
ing this feature. Hence, we will first abstract from specific 
system features such as the non-Newtonian character of 



real granular flows, non-slip boundary conditions, dissi- 
pative interactions, or force networks in quasi-static gran- 
ular flows H, 0, and discuss extensions later. This will 
allow us to show that intermittent flows are caused even 
without mechanisms like dissipative collapse ||, large 
spatio-temporal fluctuations due to force networks jfj, 
or delay- induced instabilities (as in traffic flows). These 
may magnify the effect jlfj . 

As pedestrian evacuation has been successfully de- 
scribed by driven granular particle models, where a single 
particle represents an individual pedestrian, we will for- 
mulate a common model for escaping pedestrians and 
gravity-driven outflows from vertical, two-dimensional 
hoppers. Due to the conservation of the particle number, 
we will describe the aggregate, two-dimensional particle 
flow by the continuity equation for the particle density p 
as a function of space and time. Both, the shape of a fun- 
nel and the semicircular shape of a waiting crowd suggest 
to write this equation in polar coordinates. Assuming no 
significant dependence on the polar angles 9 and ip for the 
moment, we obtain dp/dt + (l/r)d(rpv)/dr = (gener- 
alizing this to a 2d treatment later). Here, t denotes the 
time, r > the distance from the bottleneck (exit) and 
v < the velocity component in radial direction. The 
above continuity equation can be rewritten as 



dp d{pv) 
dt dr 



pv 



(1) 



where the term on the right-hand side reflects a merg- 
ing effect similar to an on-ramp or lane closure term 
in a model of freeway traffic. By use of logarith- 
mic derivatives, the above equation can be rewritten 
as dlnp(r,t)/dt = —v(r,t)dln[rp(r,t)v(r,t)]/dr. For 
the stationary case with d In p/dt = it follows from 
dhi(rpv) / dr = that the overall flow fTrrp(r)v(r) 
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through any cross section at distance r is constant: 
fTrrp(r)v(r) — firrq(r) —: — Qo — const. 



(2) 



q(r) = p(r)v(r) is the particle flow through a cross section 
of unit length. / = 1 corresponds to the half circumfer- 
ence irr of a circle of radius r, while / < 1 allows one 
to treat hoppers with an opening angle smaller than 180 
degrees. (The walls should be steeper than the angle of 
repose.) Qo > is the stationary overall particle flow. 

To facilitate the derivation of analytical results, we will 
assume the linear velocity-density relationship 



V(r) = V(p(r)) = -Vo (1 - p/Pmax) < . 



(3) 



vq means the maximum ("free") particle speed and p ma x 
the maximum particle density. Eqs. and Q give the 
quadratic equation p(r)v [l - p(r)/p max ] = Qo/(fiw) in 
p. With r crit (Qo) := <3o/(/7rg max ) it implies 



P±(r, Qo) 



1±\/1- 



:rit(Qo) 



(4) 



where (? ma x = foPmax/4 is the maximum flow. In free 
flow with d\pV (p)\j dp > 0, the density profile is deter- 
mined by the upstream boundary condition, i.e. Qo is 
given by the overall inflow Q m . Under congested con- 
ditions (d\pV(p)\/dp < 0), Qo is given by the overall 
outflow Q out = min(Q in ,2rog max ,/7rrog m ax), he. the 
minimum of the overall inflow Q; n and the maximum 
possible overall outflow. The stationary case requires 
Qin = Qout and a non- negative discriminant in Eq. 
This calls for ro > r cr i t (Qi n ), i.e. large outlets (see 
Fig.[lji). Then, p(r,t) converges to a stationary free flow 
with the density profile p_(r, Q in ) < p max /2. The den- 
sity profile for other velocity-density relationships than 
(0 can be obtained numerically. Smooth perturbations 
like the humps in Fig. QJi propagate forward at the speed 
V(p) + pdV(p)/dp = -Uq(1 - 2p/p max ) [3, compactify 
close to the outlet and leave the system. 

If the overall inflow exceeds the overall outflow (Q; n > 
Qout), particles are jammed behind the outlet (Fig. 
\3p). The initial density profile p(r, 0) again approaches 
p_(r, Qi n ) < p ma x/2 in the free-flow regime at large dis- 
tances r, but converges to p+(r, Q ou t) > Pmax/2 in the 
congested regime upstream of the outlet. The congestion 
front moves forward until the jam density p+(r, Q out ) is 
reached. Then, the particles queue up and the shock 
front at location R(t) moves backward at the speed 



dR 
~dt 



Qin Qo 



(5) 



fnR(t)[p-(R, Q in )-p+(R, Q out )] 
according to the shock wave equation 

[3 (see Fig. UJd). 

Hence, we find the free-flow density profile p(r, t) k, 
P-{r, Qin) for r > R{t), while for r < R(t) we have the 
congested density profile p(r,t) « p+(r, Q ou t)- 




r/rcrit 



FIG. 1: (Color online) Density profiles at different times, when 

(a) the inflow is low and the initial density profile has a hump, 

(b) the inflow exceeds the maximum outflow and the initial 
density profile is a step function (shock wave) . The simulation 
results have been obtained by solving the continuity equation 
with the Godunov scheme, assuming p(12r CI i t ,t) — 0.01 and 
floating boundary conditions at r — ro = 2r cr it := 2r C rit(<5in) 
in case (a), but p(2.8r crit , t) = 0.1 and Q out = 2r (?max (corre- 
sponding to the maximum outflow) with ro = 0.5r cr it in case 
(b). Note that the asymptotic density profile is p~(r, Qi n ) in 
free flow and p+(r, Q out ) in jammed flow. 



This applies to cases of continuous outflows, which 
are observed for large enough openings 0, 0] or small 
enough pedestrian velocities vq [ll|. However, if the de- 
sired velocity vq of pedestrians is high, their maximum 
density Oniax goes up and intermittent outflows are ob- 
served pi UJ|. This intermittent behavior (see Fig. EJ 
reminds of driven granular media 4] and shall be mod- 
eled now. For this, let us subdivide the particle bulk 
into shells of thickness Ar (for example, the particle di- 
ameter d or multiples of it). Within each shell of area 
A(r) w firr Ar, we assume a constant average density 
p(r,t) — N(r,t)/A(r), where N(r,t) denotes the number 
of particles in the shell of radius r at time t. Further- 
more, we assume that particles move from one shell to 
the next with velocity Vo, if they find a suitable gap, 
otherwise, they will stay. The maximum number of 
particles available to move into the shell of radius r is 
p(r + Ar, t)A(r + Ar), while the maximum number of 
available gaps in shell r is p max A(r)[l - p(r, t)/p max ], be- 
cause Pmax^4(?") is the maximum number of particles in 
the shell of radius r and q(r,t) = 1 — p(r,t)/p max rep- 
resents the fraction of free space. Finally, we assume 
that ^g(r, t) denotes the probability to find a suitable 
gap in front of a particle allowing it to move, where 
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are random numbers specified in each time step with 
< £r < 1 and ^-Ar = £r ( m or der to guarantee par- 
ticle conservation). Then, the number of inflowing par- 
ticles within the time interval At = Ar/v is N- U1 (r,t) — 
£+q(r, t) minL4(r + Ar)p(r + Ar, t), A(r)p max ] , while the 
number of outflowing particles is N out {r,t) — £r<l( r ~ 
Ar, t) mm[A(r)p(r, t),A(r — Ar)p max ]. From the balance 
equation N(r,t + At) = N(r,t)+N in (r,t)-N ont (r,t) and 
p(r,t) = N(r,t)/A(r) we get 

P{r,t) 



p{r,t + At) = P {r,t)+i+ 1- 



1 + — I p{r + Ar,t), p n 



C 1 



p(r — Ar, t) 



Pn 



pir. / ).(!-—) /» nlilx 



(6) 



Finally note that the half circle of radius ro around the 
exit is treated analogously to the shells, but we have to 
replace the area A(0) by A — nr 2 /2 and — N out (0,t) 
by — 2ropo(t)v At (i.e. the exit width 2r times the flow 
Qo{t) — Po(t)vQ) if pedestrians can leave with maximum 
velocity Vq into the uncongested space behind the exit). 
The resulting equation for the density po(t) in the last 
(sub-)area before passing the bottleneck is 



P o(t + At) = p (t)+$(t) 1- 



Po(t) 

Pmax 



x mm 



2A? 
ro 



p(r ,t),p n 



4Ar 
7rr 



Po(t). (7) 



The minimum function in Eq. 10 delineates the merging- 
related lack of space and outflow capacity. A similar 
situation and minimum function occurs in merging flows 
in urban street and production networks. With Ar = 

v At, p = P (r, t), z r = (£+ + e r -)A c = (e+ - e r -)(i - 

p/pmax)/At and for p(r, t) < (1 - Ar/r)p max we find the 
following equation in the limit At, Ar — ► 0: 



dp 



Pmax 



dp pVoCr 
dr r 



1 - 



<rP- 

(8) 

With the linear velocity-density relation this exactly 
corresponds to the previous continuity equation (Q, if 
£^ = 1, as for small enough densities (see below). Fluc- 
tuations £^ < 1, however, allow one to describe a dy- 
namics in which less particles than possible are successful 
in finding a gap in the next shell, because of coordina- 
tion problems. The random variable £^ reflects that the 
microscopic spatial configuration of the particles mat- 
ters. When the second terms in the minimum functions 
of Eq. (jSJl apply, the dynamics is given by the equation 



dp dp 

at + VoC d~r 



P{r, t)] + (rp n 



(9) 



After averaging over the noise terms , representing the 
average of £~ by £, defining the the gap density p(r, t) = 
Pmax — p{r, t), and introducing V — vq£,, this turns into a 
continuity equation for gap propagation: 



dp | djpV) 
dt dr 



pV{t) 



(10) 



Note that gaps propagate with velocity V > Q, i.e. in 
opposite direction to the particles. 

We expect that a switching between gap propagation 
and particle propagation by the minimum function can 
account for the intermittent outflows of dense granu- 
lar flows. Triggered by the randomness of the vari- 
able , the switching mechanism can produce particle 
avalanches of different sizes. The fluctuations £^ and 
their average value £ can be adjusted to experimental 
or suitable microsimulation results, e.g. to reflect the 
spatio-temporal fluctuations due to granular force net- 
works. Here, we have instead simulated Eq. © with 
binomially distributed values of t^r, i.e. = k/N with 

P(k) = ( N ) p k (l-p) N ~ k . P(k) is the probability that 



k 6 {0,1,..., N} of N ~ N(r,t) particles successfully 
manage to move forward, where p = £ is the probabil- 
ity of a particle not to be obstructed. We have used the 
phcnomcnological specification 



P(p,r) = <1 



r 

Ar 



e 7 - 



-i -n 



(11) 

if p > 0, otherwise p := 0. (/3, 7, and e are non- 
negative fit parameters.) This ensures that p(p,r) be- 
comes 1 for p — > or r — > 00 and max [0, €(7 — Ar/r)] 
for p — > pmax. That is, we have complete clogging, if 
ro/Ar < I/7 » 5/2, which reflects arching if the outlet 
2ro is too small (see Fig. EJ. Otherwise, if the density 
is low or the bottleneck far away, we have p = 1 corre- 
sponding to particles moving at the speed v = V(p). In 
queued areas with p w p max , gaps propagate upstream 
with velocity V = «o^(7 — Ar/r). 

We have presented a Shockwave approach to deter- 
mine the spatio-temporal density profile in granular 
bottleneck flows and evacuation scenarios. General- 
ized to two dimensions allowing to consider boundary 
conditions, friction, etc., there is a free-flow regime 
characterized by forward motion according to dp/dt + 
d(pV)/dx = —[pV/w(x)]dw(x)/dx, where (dw/dx)/w — 
(dV±/dx±) /V replaces 1/r and describes the relative 
change of locally available width w, i.e. the bottleneck 
effect, x is the coordinate in flow direction and V > 
the corresponding speed, while _L represents the perpen- 
dicular direction. If w(x)p > [w(x) + Axdw/dx]p max , 
the density after the next step of length Ax (where Ax 
is the mean free path or a fit parameter) would ex- 
ceed the maximum possible density p max . Therefore, if 
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Group Size 



FIG. 2: (Color online) (a) In agreement with an experiment for 
granular outflows from a two-dimensional hopper |5j, our sim- 
ulation model generates exponentially distributed avalanche 
sizes when frictionless particles with coordination problems 
are jamming at a bottleneck (i.e. theory and experiment show 
a straight line in a log-linear plot), (b) The standard devi- 
ation of the outflow, divided by the average outflow shows 
3 regimes: no outflow for ro/Ar < I/7 = 2.5, smooth out- 
flows for large outlets, and intermittent flow in between, (c) 
The relative proportion of time steps At with a stopped out- 
flow confirms this picture. Our results are quite insensitive 
to the selected parameters. For illustration, we chose — 3, 
7 = 2/5, 6 = 0.01 , Q in = 4/ At, v = Ar/At, p max = I/(Ar) 2 
and, in (a), ro = 5 Ar (jammed conditions). 

the "gap density" p(x,t) := p max — p{x,t) falls below 
the value —p max Axdw/dx, we have instead the equation 
dp/dt + d(pV)/dx = — \pV/w(x)]dw(x)/dx for gap prop- 
agation, where V < is the backward propagation speed. 
Hence, at bottlenecks we have alternating phases of for- 
ward pedestrian motion with speed V and of upstream 
dissolving pedestrian jams with average speed V, where 
V and V may fluctuate in space and time. 

These formulas are useful for the appropriate dimen- 
sioning of exits in order to avoid critical situations in 
cases of emergency evacuation of people. If the bottle- 
neck is too small (and the desired speed vo of pedes- 
trians too high), one may find intermittent flows close 
to and behind the bottleneck (Fig. b, c). These are 
due to the fact that too many "particles" are compet- 
ing for a confined space. Obviously, not all particles 
can successfully progress, when there are mutual obstruc- 
tions. This "coordination problem" has been reflected by 
a fluctuation factor the mean value p(p, w/(dw/dx)) 
of which drops significantly below 1 if p(r,t)/p max > 
1 + Ar (dw / dx) / w(x) , i.e. if not all particles fit into the 
reduced space when progressed by a distance Ax. 

The resulting dynamics is related to a stop-and-go phe- 



nomenon: In the high-density jam, the velocity is zero, 
as the particles cannot move. However, jam resolution at 
the exit causes an upstream moving shock wave, in front 
of which the density is low. Therefore, particles at the 
jam front can leave the jam. In fact, if the density in 
front of the jam is small enough, there is a forward mo- 
tion of particles filling the low-density area. Altogether, 
we will have alternating phases of jam resolution and gap 
filling processes close to the exit, which lead to alternat- 
ing propagation directions of the jam front. This may 
also explain the observed alternation in the prop agat ion 
direction of perturbations in freeway traffic flows 15]. 
Support by the DFG (He2789/7-l) is acknowledged. 
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